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Abstract 

Viral immune evasion by sequence variation is a major hindrance to HIV-1 vaccine design. To address this challenge, our 
group has developed a computational model, rooted in physics, that aims to predict the fitness landscape of HIV-1 proteins 
in order to design vaccine immunogens that lead to impaired viral fitness, thus blocking viable escape routes. Here, we 
advance the computational models to address previous limitations, and directly test model predictions against in vitro 
fitness measurements of HIV-1 strains containing multiple Gag mutations. We incorporated regularization into the model 
fitting procedure to address finite sampling. Further, we developed a model that accounts for the specific identity of mutant 
amino acids (Potts model), generalizing our previous approach (Ising model) that is unable to distinguish between different 
mutant amino acids. Gag mutation combinations (17 pairs, 1 triple and 25 single mutations within these) predicted to be 
either harmful to HIV-1 viability or fitness-neutral were introduced into HIV-1 NL4-3 by site-directed mutagenesis and 
replication capacities of these mutants were assayed in vitro. The predicted and measured fitness of the corresponding 
mutants for the original Ising model {r= —0.74, p = 3.6xl0~^) are strongly correlated, and this was further strengthened in 
the regularized Ising model {r= —0.83, p = 3.7x10"^^). Performance of the Potts model {r= —0.73, p = 9.7x10"^) was similar 
to that of the Ising model, indicating that the binary approximation is sufficient for capturing fitness effects of common 
mutants at sites of low amino acid diversity. However, we show that the Potts model is expected to improve predictive 
power for more variable proteins. Overall, our results support the ability of the computational models to robustly predict 
the relative fitness of mutant viral strains, and indicate the potential value of this approach for understanding viral immune 
evasion, and harnessing this knowledge for immunogen design. 
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Introduction 

The ideal way to combat the spread of HIV-1 is with an 
effective prophylactic or therapeutic vaccine [1,2]. One of the 
greatest challenges hindering the achievement of this goal is the 
incredible sequence diversity and mutability of HIV- 1 [3] , which 
can limit the effectiveness of the immune response [2,4]. 

CD8-I- T cells are instrumental in reducing viral load in HIV-1 
acute infection [5] and in maintaining the viral set point during 
chronic HIV/SIV infection [6,7]. However, HIV-1 is able to 
escape the CD8-I- T cell response through mutations in or adjacent 
to HIV- 1 epitopes that are presented by HLA class I molecules on 



the surface of the infected cells [7]. One proposed strategy for 
reahzing a potent prophylactic or therapeutic vaccine is to target 
CD8-I- T cell responses to conserved regions of HIV-1, aiming to 
reduce incidences of immune escape or, if escape occurs, to reduce 
viral fitness and lower the viral set point, thereby slowing disease 
course and reducing transmission at the population level [8,9]. 
While escape mutations at highly conserved sites often damage the 
viability of virus [10], this approach is confounded by the 
development of compensatory mutations which restore or partially 
restore viral fitness [9] . Thus, to maximize the effectiveness of a 
vaccine-induced immune response one must look beyond conser- 
vation of single residues to identify regions where mutations are 



PLOS Computational Biology | www.ploscompbiol.org 



1 



August 2014 | Volume 10 | Issue 8 | el 003776 



The Fitness Landscape of HIV-1 Gag 



Author Summary 

At least 70 million people have been infected with HIV 
since the beginning of the epidemic and an effective 
vaccine remains elusive. The high mutation rate and 
diversity of HIV strains enables the virus to effectively 
evade host immune responses, presenting a significant 
challenge for HIV vaccine design. We have developed an 
approach to translate clinical databases of HIV sequences 
into mathematical models quantifying the capacity of the 
virus to replicate as a function of mutations within its 
genome. We have previously shown how such "fitness 
landscapes" can be used to guide the design of vaccines to 
attack vulnerable regions from which it is difficult for the 
virus to escape by mutation. Here, using new modeling 
approaches, we have improved on our previous models of 
HIV fitness landscape by accounting for undersampling of 
HIV sequences and the specific identity of mutant amino 
acids. We experimentally tested the accuracy of the 
improved models to predict the fitness of HIV with 
multiple mutations in the Gag protein. The experimental 
data are in strong agreement with model predictions, 
supporting the value of these models as a novel approach 
for determining mutational vulnerabilities of HIV-1, which, 
in turn, can inform vaccine design. 



not only highly deleterious, but where further mutations elsewhere 
in the proteome are unlikely to restore lost fitness, but rather, lead 
to additional fitness costs due to deleterious synergistic efiects. 

Our group has developed computational models to identify such 
vulnerable regions of the HIV-1 proteome and to predict the 
fitness landscape of HIV- 1 proteins, providing tools for designing 
vaccine immunogens that may limit both HIV- 1 evasion of CD8-H 
T cell responses and the development of compensatory mutations 
[1 1,12]. In an early qualitative study we identified groups of amino 
acids in HIV-1 Gag coupled by structural and functional 
constraints that cause these residues to co-evolve with each other, 
but evolve nearly independently of the other residues in the 
protein [12]. In analogy with past studies on the economic markets 
and enzymes [13-15], we termed these groups of residues 
"sectors". This analysis and human clinical data revealed one 
sector in Gag, which we termed sector 3, where multiple mutations 
were more likely to be deleterious. This group of residues is 
naturally targeted more by elite controllers [12]. It is expected to 
be particularly vulnerable to CD8-H T cell responses that target 
multiple residues in it since multiple mutations within this sector 
are likely to significandy diminish viral fitness, thereby restricting 
available escape and compensatory paths [12]. 

This approach, however, does not allow us to determine 
precisely which residues should be targeted, as it does not quantify 
the relative replicative viability of viral strains bearing specific 
mutations. Nor does it identify viable escape routes that remain 
upon targeting residues in the vulnerable regions, or inform how 
best to block them. To begin to address these issues, we developed 
a computational model, rooted in statistical physics, which aims to 
predict the viral fitness landscape (viral fitness as a function of 
amino acid sequence) from sequence data alone and applied it to 
HIV-1 Gag [1 1]. Similar methods have previously been employed 
to study other complex biological systems, from describing the 
activity patterns of neuronal networks [16-19] to the prediction of 
contact residues in protein families [17,20,21]. 

The idea underlying our approach is to first characterize the 
distribution of sequences in the population, which we expect to be 
correlated with fitness (see below). Due to the small number of 



available sequences compared to the size of the sequence space, 
direct estimation of the probability distribution characterizing the 
available sequences is precluded. Thus, we instead aim to infer the 
least biased probability distribution of sequences that fits the 
observed frequency of mutations at each site, and all correlations 
between pairs of mutations (the one- and two-point mutational 
probabilities). Mathematically, "least biased" implies the distribu- 
tion that has maximum entropy in the information-theoretic sense 
[22]. The maximum entropy distribution that fits the one- and 
two-point mutational probabilities has a form reminiscent of that 
describing equilibrium configurations of an Ising model in 
statistical mechanics. We generated such models using multiple 
sequence alignments (MSA) for the four subunit proteins of Gag in 
HIV-1 clade B [1 1] (described in Supporting Information Text SI, 
Section 1). This model assigns to each viral strain an "energy" (E), 
which is inversely related to the probability of observing this 
sequence. 

We expect more prevalent sequences to be more fit, consistent 
with expectations from simple models of evolution [23] though the 
precise correspondence between fitness and prevalence may have 
a more compUcated dependence on factors such as the shape of 
the fitness landscape, as predicted by quasispecies theory [24]. 
Furthermore, this expectation cotdd be confounded by immune 
responses in the patients from whom the virus samples were 
collected, and phylogeny. Recent analyses suggest (described more 
fully in the discussion) that in spite of these effects, at least for Gag 
proteins, the rank order of prevalence and in vitro replicative 
fitness should be similar [25]. Strains with high E values are 
predicted to be less fit tiian strains with low E values. Predictions of 
the model seemed to be in good agreement with experimental data 
on in vitro replicative fitness, as well as clinical observations on the 
frequency and impact of viral escape mutations [11]. 

Our aim in the current work is twofold. First, we present new 
advances in the inference and modeling of viral fitness landscapes 
that address previous theoretical and computational limitations. 
Second, we describe new in vitro fitness measurements for viruses 
containing multiple Gag mutations, performed to further test 
fitness predictions using the improved computational methods. To 
give a broad test of the predictive power of the fitness models, we 
have performed comparisons for HIV-1 strains containing 
multiple mutations predicted to harm HIV-1 viability as well as 
combinations predicted to be relatively fitness neutral. We find 
that fitness measurements of these mutant strains are in good 
agreement with model predictions. 

Methods 

Computational models to translate sequence data to 
viral fitness landscapes 

Our key hypothesis in formulating models of HIV fitness is that 
the prevalence of viruses with a given sequence, that is, how often 
the sequence is observed, is related to its fitness. Simply, fitter 
viruses should be more frequent in the population than those that 
are unfit. This hypothesis can be proven for some idealized 
evolutionary models [23], but cannot be made exact for the 
complicated nonequilibrium host-pathogen riposte between hu- 
mans and HIV. However, our theoretical work, backed by 
extensive computational studies, suggests that the rank order of 
fitness and prevalence of strains should be strongly monotonicaUy 
correlated, provided we compare sequences that are phylogenet- 
ically close [25]. Thus, if we construct a model to predict the 
likelihood of observing different viral strains with given sequences, 
it can predict the relative fitness of the strains. We achieved this 
goal by constructing a maximum entropy model for the 
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probability of observing sequences in the MSA [26] . The simplest 
model in this class is an Ising model, a simple model of interacting 
binary variables from statistical physics which has been widely 
applied to study collective behavior in complex systems. The 
parameters of this Ising model are obtained by imposing the 
constraint that it reproduce the pattern of correlated mutations 
(relative to the consensus sequence) observed in a multiple 
sequence alignment (MSA) of HIV-1 amino acid sequences 
extracted from infected hosts. Specifically, the parameters were 
chosen such that the frequency of mutations at each single residue 
and the frccjucnc)" of simultaneous mutations at each pair of 
residues were the same in both the Ising model and the MSA. 
Importantly, the model also reproduced higher order mutational 
correlations accurately, even though these mutational frequencies 
were not directly fitted [11]. 

As described in our previous publication [1 1], in the Ising model 
amino acid sequences in the MSA are compressed into binary- 
strings by assigning a 0 to each position where the amino acid 
matches the consensus sequence ("wild-type"), and a 1 to each 
position with a mismatch ("mutant"). While this binary approx- 
imation greatly simplified our modefing approach, the reduction in 
complexity has several drawbacks. Firstly, there is a loss of residue- 
spc'cific resolution. The fitn("ss predictions of our model are 
insensitive to tlie precise identity of mutant amino acids, and thus 
the model cannot resolve fitness differences between proteins 
containing different mutant amino acid residues in a particular 
position. Secondly, for relatively conserved proteins such as HIV-l 
Gag, where the number of viable amino acids at each position is 
rather limited, this binary simplification represents a reasonable 
approximation. However, it is less justified for highly mutable 
proteins where the wild-type residue in each position is not the 
overwhelmingly most probable amino acid, as is the case for the 
HIV-1 Env protein. 

In our original approach, we fit the Ising model parameters to 
precisely reproduce the observed one and two-residue mutational 
correlations within the MSA. However, simultaneous mutations at 
certain pairs of residues were never observed. This led to another 
deficiency in our original modeling approach in that pairs of 
mutations not observed in the MSA were predicted to be completely 
unviable (E = oo). While it is possible that such mutant viral strains 
have exacdy zero replicative fitness, it is more likely that they are 
highly unfit strains (possessing non-zero replicative fitness) that 
simply arise too seldom to be observed within our finite-sized MSA. 

In this work, we present three significant advances of our 
original model to predict viral fitness, which also the aforemen- 
tioned limitations. First, we incorporate Bayesian regularization 
into our fitting procedure to eliminate the prediction of zero 
replicative fitnesses for mutations not present within our MSA. 
Second, we implement a new algorithm for inferring an Ising 
model from sequence data, which dramatically accelerates the 
computation of model parameters. Third, we relax the binary 
approximation to infi'r \'iral fitness landscapes that explicitly retain 
the amino acid identities at eac;h position. We achieve this by 
describing the viral fitness landscape using a multistate general- 
ization of the Ising model known as the Potts model, another 
established and well-studied model in statistical physics [27]. We 
also implement Bayesian regularization into the fitting of the Potts 
model parameters. 

Model 1: Regularized and computationally fast inference 
of Ising models of viral fitness 

Inference of the parameters of the Ising models, commonly 
referred to as the inverse Ising problem, is a canonical inverse 
problem lacking an analytical solution that may be tackled in 



many ways [16,17,19-21,28,29]. We improve upon our previous 
techniques described in [1 1] by incorporating regularization and 
implementing new inference algorithms, which gready decrease 
the computational burden and accelerate model fitting. 

To control the effects of undersampling and to improve the 
predictive power of the inferred fitness models, we incorporate 
Bayesian regularization into our inference algorithm [18,19,30,31] 
in the form of a Gaussian prior distribution for the model 
parameters describing pairwise couplings between residues (see 
Text SI, Sections 1.3 and 2.5). Regularization of this form is also 
known as Tikhonov regularization or ridge regression [32]. With 
this addition, the probability of observing any sequence, including 
those containing pairs of mutations not observed in the MSA, is 
nonzero. We have also computed a correction to the energ\r of 
each sequence to account for the possible bias that strains near 
fitness peaks are more likely to be observed than would be 
expected from their intrinsic fitness when sampled from a finite 
distribution (see Text SI, Section 3.2). 

In an algorithmic advance over our previous fitting procedure, we 
fit the parameters of our regularized Ising model using the selective 
cluster expansion algorithm of Cocco and Monasson [18,19] which 
identifies clusters of strongly interacting sites and iteratively builds a 
solution for the whole system by solving the inverse Ising problem 
for each cluster. With this approach, we cut the CPU time necessary 
to infer the parameters of the Ising model from roughly 1 2 years 
[1 1] to 5 hours for p24, an improvement by four orders of 
magnitude. Roughly, we expect algorithm run-time to scale as 0{Nn 
exp(n)), where is the system size (number of amino acids) and n is 
the size of a typical "neighborhood" of strongly interacting sites. For 
a review and appfications of this method see [18,30]. Complete 
details of our modeling approach and numerical fitting procedures 
are provided in the Text S 1 , Section 1 . 

Model 2: Regularized Potts models of viral fitness 

An ideal model of viral fitness would be able to capture the fuU 
(unknown) distribution of correlated mutations throughout the 
sequence, and thus reproduce the prevalence of every viral strain. 
Sequences in the MSA represent a sample of the possible strains of 
the virus, providing information about the distribution of point 
mutations, pairs of simultaneous mutations, triplets of simulta- 
neous mutations, and all higher orders. However, since the 
number of available sequences in the MSA is very small compared 
to the size of the accessible sequence space, and because mutations 
at most sites are rare, higher order mutations will be severely 
undersampled. Thus, following our previous approach we appeal 
to the maximum entropy principle to seek the simplest possible 
model capable of reproducing the single site and pair amino acid 
frequencies [11,22], for which the problem of undersampling is 
less severe. From this analysis, the Potts model is the least 
structured model capable of reproducing the one and two-position 
frequencies of amino acids observed within the MSA [21]. 

To introduce the Potts model, we represent the sequence of a 
particular m-residue protein as a vector, A = {Ai^Yk=\7 where the 
elements Ai^ can take on the q — 2l integer values [1, 2, 21] 
denoting an arbitrary encoding of the 20 natural amino acids, plus 
a gap [21]. In the Potts model the probability PiA) of observing a 
particular sequence A is given by 
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In analogy with the statistical physics literature, we refer to E as 
a dimensionless "energy," the function E{A) as the Hamiltonian, 
and the normalizing factor Z as the partition function [27]. The 
model is paramc't(;rized by a set of m q-dimensional vectors, 
{hi}'^^i, and a set of m(OT— 1)/2 q-by-q matrices, {•/i/},<^-- The hi 
vectors give the contribution of the identity of each amino acid in 
each position to the overall sequence energy, and the fy matrices 
give the contribution to the energy of pairwise interactions 
between amino acids in different positions. 

To fit the Potts model, we implemented a generalization of the 
semi-analytical extension of the iterative gradient descent imple- 
mented by Mora and Bialek [11,17]. This approach implements a 
multi-dimensional Newton search to iteratively adjust the 
model parameters until the predictions of the model for the one 
and two-position frequencies of amino acids reproduce those 
observed within the MSA. In an advance over the original 
incarnation of this algorithm, we have derived closed form 
expressions for the gradients required by the Newton search, 
thcrcb)' ol)viating the need for their numerical estimation by finite 
differences (which would result in a more computationally 
expensive and less numerically stable secant search procedure). 
Our approach is semi-analytical in the sense that while we have 
analytical expressions for the Newton search gradients, we use a 
Monte Carlo procedure to numerically estimate the one and two- 
position amino acid frequencies predicted by the model at each 
stage of parameter refinement. We are currentiy developing a 
Potts generalization of the cluster expansion algorithm [19] to 
accelerate fitting. We incorporate Bayesian regularization into our 
fitting procedure in a precisely analogous manner to that described 
above for the regularized Ising model by introducing a Gaussian 
prior distribution over the J,j parameters. Inference of the Potts 
model parameters for p24 required approximately 1.4 years of 
CPU time using a generalization of the gradient descent approach 
described in Ref [11]. Fitting the model parameters by gradient 
descent is expected to scale as ()[[Np'f), where is the number of 
amino acids in the protein, and /; is the characteristic number of 
mutant residues observed at each position. Full details of the fitting 
and regularization procedures are provided in Text SI, Section 2. 
The code implementing the inverse Potts inference algorithm is 
also provided in Supporting Information Code SI. 

In vitro experiments 

To test the accuracy of these models in predicting the fitness 
landscape of HIV-1 Gag, we performed in vitro experiments to 
measure the fitness of various Gag mutants. Previously we had 
measured the in vilro replication capacities of 19 Gag p24 
mutants, 16 of which contained single mutations in Gag p24, and 
compared these with fitness predictions of our original Ising model 
[11]. Here, we extend this work to measure the replication 
capacities of HIV-1 strains containing various combinations of 
mutations, predicted to be either harmful to HIV-1 viability or 
fitness-neutral, in Gag p24 and pi 7 and we compare measure- 
ments not only to the original Ising model described in Ref [1 1], 
but also to regularized versions of Ising and Potts models that we 
have developed here. Specifically we considered 17 mutations 
pairs, one triple, and 25 single mutations within these combina- 
tions, as fisted in Table 1. These mutations were introduced into 
the widely used laboratory-adapted HIV-1 clade B reference strain 
NL4-3. 

The tested mutants can be divided into 4 categories, viz. (i) Gag 
p24 pairs with high E values located within a group of co-evolving 
amino acids termed sector 3 (cf Ref [12]), (ii) HLA-associated Gag 
p24 pairs with high E values, (iii) Gag p24 pairs/ triple with low E 



values, and (iv) Gag pi 7 pairs (Table 1). These mutation 
combinations were chosen according to E values predicted by 
the published Ising model [11], where -E>90 or £ = oo were 
considered high E values and £< 15 were considered low E values. 
Note that, due to the couplings between mutations at different 
sites, parameterized by the Jy in equation 1 , the E values depend 
not only on the specific mutations introduced but also on the 
sequence background. The E values for mutations reported here 
are computed with the HIV-1 NL4-3 sequence background, which 
differs from the pi 7 and p24 MSA consensus sequences by 8 
mutations (R15K, K28Q, R30K, K76R, V82I, T84V, E93D, 
S125N1 and 2 mutations (N252H, A340G), respectively. The p24 
region of Gag was focused on since this is the most conserved 
region of the protein. First, we selected six mutation pairs, 
predicted to be unfavorable in combination, in sector 3 of Gag p24 
since we previously found this to be an immunologically 
vulnerable group of co-evolving residues in which multiple 
mutations are not well-tolerated [12]. Since it is desirable to 
identify low fitness/ non-viable combinations of escape mutations 
for vaccine immunogen design aimed at reducing viral fitness or 
blocking viable escape pathways, we aimed to identify pairs of 
likely escape mutations with high E values. Virus mutations that 
are statistically associated with the expression of specific host HLA 
class I alleles, which also restrict the same epitopes in which the 
mutations are found, are likely to be CDB-I- T cell-driven escape 
mutations [33]. We therefore tested five high E pairs of mutations 
located at HLA-associated Gag p24 codons (HLA-associated 
variants defined in [34,35]) in or next to optimal CD8+ T cell 
epitopes (A-list epitopes from the Los Alamos HIV sequence 
database [36]) that were restricted by the same HLA. For 
comparison with high E mutation pairs, mutation combinations 
with low predicted E values were included in testing, comprising 
known favorable compensatory pairs in Gag p24 where 219Q^ 
compensates for the 242N escape mutant [37] and 147L 
compensates for the 146P escape mutant [10], as well as one 
pair in sector 3 of Gag p24 and a Gag p24 triple mutant. 
Additionally, for broader testing, two mutation pairs in Gag pi 7 
were selected. We note that the most commonly observed mutant 
amino acid at each codon was tested. 

We introduced these mutation combinations into the HIV-1 
NL4-3 plasmid by site-directed mutagenesis and their presence was 
confirmed by sequencing, as described previously [38]. Generation 
of mutant viruses from mutated plasmids and the measurement of 
their replication capacities were performed as previously [11,38]. 
Briefly, mutated plasmids were electroporated into an HIV-1- 
inducible green fluorescent protein reporter T cell line, harvested at 
—30% infection of cells, and the replication capacities of the 
resulting mutant viruses were assayed by flow cytometry using the 
same ceU line. Replication capacities were calculated as the 
exponential slope of increase in percentage infected cells from days 
3-6 following infection at a MOI of 0.003, normalized to the growth 
of wild-type NL4-3 (i?C=l). Three independent measurements 
were taken and averaged. Mutant viruses were re-sequenced to 
confirm the presence of introduced mutations. 

Results 

Comparison of predictions from different models 

The values of E predicted by our original and new modeling 
approaches for the 43 HIV-1 NL4-3 Gag mutants tested here are 
shown in Table 2. Absolute comparison of the E values between 
the models are not meaningful, but the relative E values of mutants 
are generally in excellent concordance between models (Pearson's 
correlation, r>0.85 and^S5.3xlO~", two-tailed test). 
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Table 1. HIV-1 NL4-3 Gag mutants selected for testing of predicted energy costs (E) by an in vitro HIV-1 replication capacity assay. 



Mutant 


Gag subunit 


Category of pairs/triple 


E" 


1861 


p24 




78.74 


269E 


p24 




43.43 


1 86I269E 


p24 


Sector 3", high E" 


Infinity 


295E 


p24 




22.81 


186I295E 


p24 


Sector 3, high E 


Infinity 


181R 


p24 




44.62 


310T 


p24 




6.26 


181R310T 


p24 


Sector 3, high E 


Infinity 


182S 


p24 




25.13 


198V 


p24 




Infinity 


182S198V 


p24 


Sector 3, high E 


Infinity 


179G 


p24 




56.09 


229K 


p24 




44.63 


179G229K 


p24 


Sector 3, high E 


97.01 


174G 


p24 




Infinity 


243 P 


p24 




66.65 


1 74G243P 


p24 


Sector 3, high E 


Infinity 


1681 


p24 




38.58 


31 5G 


p24 




19.11 


168I315G 


p24 


HLA-associated, high E 


Infinity 


331R 


p24 




11.77 


1 861331 R 


p24 


HLA-associated, high E 


Infinity 


302R 


p24 




11.10 


302R315G 


p24 


HLA-associated, high E 


Infinity 


315G331R 


p24 


HLA-associated, high E 


Infinity 


1901 


p24 




41.52 


190I302R 


p24 


HLA-associated, high E 


Infinity 


21 9Q 


p24 




6.73 


242N 


p24 




8.68 


219Q242N 


p24 


p24, low E, compensatory 


10.80 


146P 


p24 




7.22 


147L 


p24 




3.42 


146P147L 


p24 


p24, low E, compensatory 


6.58 


326S 


p24 




4.59 


310T326S 


p24 


p24, low E, sector 3 


10.53 


173T 


p24 




5.92 


1 73T286K 


p24 




4.89 


173T286K147L 


p24 


p24, low E, triple 


4.12 


12K 


pi 7 




3.74 


12K54A 


pi 7 


pi 7, low E 


4.84 


86F 


pi 7 




8.00 


92M 


pi 7 




8.74 


86F92M 


pl7 


pi 7, high E 


Infinity 



^Energy cost predicted by original Ising model [1 1] taking into account differences in NL4-3 and the multiple sequence alignment used in model generation. E is 2.98 for 
wild-type NL4-3 p24 and 3.43 for wild-type NL4-3 pi 7. 

''Mutation pairs within an immunologically vulnerable group of co-evolving residues, termed sector 3, that we previously identified qualitatively [12]. 
'^f>90 or were considered high E values and E<15 were considered low E values. 
doi:l 0.1 371 /journal.pcbi.1 003776.t001 
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Table 2. Energy costs (E) of HIV-1 NL4-3 Gag mutants predicted by computational models. 




Mutant 


Gag subunit Ising E* 


Regularized Ising E** 


Regularized Potts E*^ 


1861 


p24 


78.74 


9.98 


11.24 


269E 


p24 


43.43 


11.46 


12.18 


186I269E 


p24 


Infinity 


17.77 


18.97 


295E 


p24 


22.81 


9.05 


11.03 


186I295E 


p24 


Infinity 


15.36 


17.79 


181R 


p24 


44.62 


13.55 


12.12 


31 OT 


p24 


6.26 


5.87 


7.20 


181R310T 


p24 


Infinity 


15.74 


14.87 


182S 


p24 


25.13 


7.11 


9.68 


198V 


p24 


Infinity 


12.32 


_d 


182S198V 


p24 


Infinity 


15.77 


_d 


179G 


p24 


56.09 


11.14 


11.57 


229K 


p24 


44.63 


10.52 


11.68 


179G229K 


p24 


97.01 


17.99 


18.81 


174G 


p24 


Infinity 


15.47 


11.71 


243P 


p24 


66.65 


11.1 


11.08 


1 74G243P 


p24 


Infinity 


22.9 


18.32 


1681 


p24 


38.58 


9.8 


10.30 


31 5G 


p24 


19.11 


6.85 


10.64 


168I315G 


p24 


Infinity 


14.78 


16.39 


331R 


p24 


11.77 


7.37 


9.17 


1 861331 R 


p24 


Infinity 


13.68 


15.85 


302R 


p24 


11.1 


7.75 


9.23 


302R315G 


p24 


Infinity 


12.4 


15.30 


315G331R 


p24 


Infinity 


10.56 


15.22 


1901 


p24 


41.52 


8.2 


11.41 


190I302R 


p24 


Infinity 


12.28 


16.12 


21 9Q 


p24 


6.73 


5.65 


6.90 


242N 


p24 


8.68 


6.7 


8.05 


219Q242N 


p24 


10.8 


8.04 


10.07 


146P 


p24 


7.22 


5.62 


6.26 


147L 


p24 


3.42 


4.25 


6.54 


146P147L 


p24 


6.58 


5.77 


4.74 


326S 


p24 


4.59 


4.78 


5.69 


310T326S 


p24 


10.53 


7.72 


8.81 


173T 


p24 


5.92 


5.81 


7.02 


1 73T286K 


p24 


4.89 


6.56 


7.75 


173T286K147L 


p24 


4.12 


5.93 


6.78 


12K 


pi 7 


3.74 


1.91 


4.38 


12K54A 


pl7 


4.84 


3.19 


5.63 


86F 


pi 7 


8 


4.53 


6.00 


92M 


pl7 


8.74 


6.01 


9.43 


86F92M 


pl7 


Infinity 


9.52 


12.57 


"E is 2.98 for wild-type NL4-3 p24 and 3.43 for wild-type NL4-3 pl7. 
""E is 3.67 for wild-type NL4-3 p24 and 1.64 for wild-type NL4-3 pi 7. 
'E is 4.43 for wild-type NL4-3 p24 and 2.81 for wild-type NL4-3 pi 7. 

^The 198V mutation was not observed within the MSA used to fit the Potts model, precluding the fitted model from assigning an energy to viral strains containing this 
point mutation. 

doi:l 0.1 371/journal.pcbi.l 003776.t002 
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Figure 1. Replication capacities of NL4-3 viruses encoding mutations in HIV-1 Gag. Graphs show replication capacities of NL4-3 viruses 
encoding (A) Gag p24 mutation pairs with high E values that were previously identified to be in vulnerable co-evolving groups [12] and single 
mutations within these pairs; (B) Gag p24 HLA-associated pairs with high E values and single mutations within these pairs; (C) Gag p24 pairs/triple 
with low E values as well as single mutations within these combinations; and (D) Gag pi 7 pairs including single mutations within the pairs. Those 
mutants that (i) were not viable or (ii) were not viable unless further mutations developed (indicated with an asterisk), were assigned a replication 
capacity of zero. Mutation pairs and triples are shown in grey while single mutations within these combinations are shown in black. Replication 
capacities of mutant viruses are expressed relative to the replication capacity of wild-type NL4-3 virus (/?C= 1). Bars represent the mean of three 
independent experiments and error bars represent standard deviation from the mean. 
doi:1 0.1 371 /journal.pcbi.1 003776.g001 



Experimental findings 

The in vitro fitness measurements for all mutants, grouped 
according to categories, are shown in Figure 1. We initially 
compared our model predictions and fitness measurements for 
each category of mutant pairs to evaluate whether mutant 
combinations with high and low predicted E values corresponded 
to substantial fitness cost or litde/ no fitness cost, respectively. 

Briefly, all Gag p24 sector 3 mutation pairs with high E values 
were not viable in our assay system, and were assigned a 
replication capacity of zero (Figure lA). Similarly, with the 
exception of 315G331R, the ii\ c- high E HLA-ass<K:iated mutation 
pairs showed substantial ri^duction in replication capacity, to 
between 0-56% of wild-type levels (Figure IB). Non-viable 
mutants (RC = 0) were those for which the generation of virus 
stocks from plasmids encoding these mutation pairs failed, or, in 
two instances - mutants 186I295E and 1 86133 IR - were not 
viable unless further mutations developed, confirming unfavor- 
abUity of the mutation combination. Briefly, concentrated virus 
stocks for mutants 186I295E and 1 86133 IR were harvested at > 
22 days post-electroporation compared with the median har\'est- 
ing time of 6 days post-electroporation for all mutants (at which 
time the 186I295E and 186I331R mutants had infected =1% 
cells). Sequencing of these viruses revealed the presence of 
additional mutations and/ or reversion of introduced mutations. 
For mutant 186I295E, amino acid mixtures were detected at 
codons 63 (Q/R), 177 (D/E) and 186 (I/V), and for mutant 
186133 IR, mixtures were detected at codons 168 (I/V) and 331 
(K/R), as well as reversion of 1861 to 186T. On repeating virus 
generation for these mutants, additional mutations similarly 
developed - mixtures were observed at codons 214 (R/K) and 
271 (N/S) for mutant 186I295E, and 232 (R/M) and 260 (D/E) 
for mutant 186I331R. With the exception of 1861295E and 
186133 IR, sequencing confirmed that all mutant viruses had only 
the specific mutations introduced. The spontaneous mutations 
186V, 271S and 232M were not observed in die MSA and die 
new mutation combinations did not have lowered E values in any 
of the models, with the exception of tlu- incomplete 
186133 1R260D combination (complete observed combination 
1861, 33 IR, 232R/M, 260D/E) which displayed a slighdy lower 
energy than 1 86133 IR in the regularized Ising model only (11.5 
vs. 13.7) (data not shown). Nevertheless, these observations 
confirm that 186I295E and 186I331R are unfit mutation 
combinations requiring compensatory paths to restore viability. 
Taken together, the data on high E p24 mutants confirm mutation 
combinations predicted to be unfit, and also identify combinations 
of HLA-associated mutations in/next to optimal CD8-I- T cell 
epitopes (mutations likely to result in CD8-I- T cell escape [33]) that 
carry substantial fitness costs. 

Those p24 mutation combinations, including known compen- 
satory pairs, that were predicted to have low E values displayed 
replication capacities similar to that of wild-type NL4-3, indicating 
that these combinations had littie or no cost to HIV- 1 replication 
capacity in accordance with predictions (Figure IC). Similarly, all 
pi 7 mutants tested had replication capacities close to that of the 



wild-type NL4-3 virus, consistent with the predicted E values of all 
mutants except 86F92M (Figure ID). 

Overall, for only two (86F92M and 315G331R) of die 17 
mutant pairs the fitness measurement did not correspond to the E 
value prediction of high or low fitness cost. It should however be 
noted that the disparity between E values and measured 
replication capacities for these mutant pairs is somewhat mitigated 
in the regularized models. The E values for the regularized Ising 
model for these mutants (which were assigned an E value of 
infinity by the original Ising model) are lower than those of other 
mutants previously assigned infinite energies, and the same is true 
for mutant 86F92M in the regularized Potts model. 

Quantitative comparison between in silico predictions 

and in vitro measurements 

Next, we assessed the relationship between fitness measure- 
ments and E values predicted by our original Ising, regularized 
Ising and regularized Potts models using Pearson's correlation 
tests. There is a strong correlation between the metric of fitness 
(values of E, Table 1) predicted by the original unregularized Ising 
model and our experimental measurements (Pearson's correlation, 
r = -0.74 andp = 3.6 xlO"**, two-tailed) (Figure 2A), however this 
correlation out of necessity excludes mutants with E values equal 
to infinity (n=13). The regularized Ising model allows for 
inclusion of these data points resulting in a stronger correlation 
between predictions and fitness measurements (Pearson's correla- 
tion, r= -0.83 and p = 3.7x10"'^, two-taHed) (Figure 2B), which 
is slighdy improved by focusing on Gag p24 mutants only 
(Pearson's correlation, r = —0.85 a.ndp = 1.4xl0~", two-tailed). 
There is also a strong agreement between the residue-specific Potts 
model energies and replication capacity (Pearson's correlation, 
r = -0.73 and p = 9Jx 10"'', two-tailed) (Figure 2C). 

In practice, one may be concerned with a more coarse-grained 
measure of viral fitness: will a virus with a given sequence be able 
to replicate with similar efficiency to the wild-type, or will it be 
significantiy impaired? To explore this point, we grouped the 
experimentally tested mutants into two categories, "fit" {RC^O.5) 
and unfit (RC<0.5), and tested the ability of the fitness landscape 
models to predict which class each sequence would belong to 
based on their E values. This was accomphshed by fitting a hnear 
classifier to the data using logistic regression (Text S 1 , Section 3. 1). 
The regularized Ising model E classifier is highly accurate (9 1 % 
accuracy at optimal threshold, AUROC = 0.93) - we observed a 
strong, significant difference in replication capacities between the 
mutants classified as unfit and those classified as fit (Mann- 
Whitney U=32, p = 4.5xl0-^) (Figure 3 A). Specifically, four 
mutants (86F92M, 1901, 190I302R and 243P) were not classified 
correctiy. However, 190I302R, which was classified as unfit 
{E = 8.6), exhibited a fitness close to that of the 0.5 cutoff 
{RC = 0.56) and 243P, which displayed low fitness (iJC = 0.36), 
had a predicted E value {E = 7.4) bordering on the classifier E 
value. The Potts model classifier also performs well (81% accuracy 
at optimal threshold, AUROC = 0.80), but provides a slighdy 
weaker difference between the fit and unfit classes (Mann- Whitney 
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Figure 2. Relationship between predicted E values and 
replicative capacities of HIV-1 NL4-3 Gag mutants. Scatter plots 
showing strong correlations between measured replication capacities 
of mutants and E values predicted by (A) original Ising (Pearson's 
correlation, /•=-0.74 and /7 = 3.6x 10"'', two-tailed test, n = 30), (B) 
regularized Ising (Pearson's correlation, r= —0.83 and p = 3Jx 10"'^, 
two-tailed test, n = 43) and (C) regularized Potts (Pearson's correlation, 
r=-0.73 and /) = 9.7xl0"', two-tailed test, n = 41) models. In the 
original Ising model (panel A), mutants with E values of infinity (n = 13) 
are excluded from the correlation. 
doi:1 0.1 371/journal.pcbi.l 003776.g002 

U = 70, p = 4.0 X 10^^) (Figure 3B). Here, seven mutants were not 
classified correctly, including the same four not classified correctly 
by the regularized Ising model as well as mutants 174G, 181R, 
269E and 315G331R. Similar to mutant 243P, mutants 174G, 



181R and 269E were unfit {RC = 0) but had a predicted E values 
ranging from 7.3 to 7.7, fairly close to that of the classifier E value. 

Discussion 

In this study, we have substantially advanced our modeling 
approaches and tested the predictive power of these models by in 
vitro fitness measurements of HIV encoding various mutation 
combinations in the Gag protein. The in vitro functional data are 
overall in strong agreement with the viral fitness landscape models 
and support the capacity of these models to robustly predict both 
continuous and "coarse-grained" measures of HIV-1 in vitro 
replicative fitness. Performance of the regularized Potts and 
regularized Ising models here is similar, which is not unexpected as 
Gag in general is not highly mutable and the mutants tested here 
were the most common ones, making the binary approximation a 
fairly good assumption. Indeed, in instances where the binary 
approximation is valid, we might encounter poorer performance 
from the Potts model relative to the Ising due to a diminished ratio 
of samples (i.e., sequences in the MSA) to parameters (i.e., h and / 
values) making robust numerical fitting of the former more 
challenging than the latter. It is nevertheless encouraging that we 
are capable of fitting a significantly more complicated Potts model 
that retains residue-specific resolution without compromising the 
fidelity of our predictions. Improved inverse Potts inference 
methods which better meet these numerical challenges may also 
improve performance of the Potts model with respect to the Ising 
model results. 

Simple theoretical analysis suggests that models which differen- 
tiate between different mutant amino acids at the same site, like 
the Potts model employed here, wiU be necessary to make fitness 
predictions for highly mutable proteins such as Env and Nef, or to 
predict the fitness of sequences containing sites with mutations to 
less frequendy observed amino acids. Using a simple toy model, we 
show in Text SI, Section 3.4 that the binary approximation (Ising 
model) has several potential deficiencies compared to a Potts 
model. In particular, the Ising model generically overestimates the 
fitness of mutant sequences, particularly for sequences containing 
uncommon mutations. Also, in the Ising case the inferred 
interaction between mutations at different sites is dominated by 
the interaction between the most common mutants, while the Potts 
model is able to accurately capture interactions between rare 
mutants. Future work wiU involve testing Ising and Potts model 
predictions for more highly variable proteins and for mutations to 
uncommon amino acids. 

While this study confirms the usefulness of this method for 
predicting HIV-l replicative fitness, at least for closely related 
sequences, caution wiU be necessary in applying this method to 
predict the relative fitness of multiple strains separated by a large 
number of mutations. In the measure of prevalence used to infer 
the Ising and Potts model fitness landscapes, factors such as 
phylogeny are implicitly included. Analysis conducted in [25] 
suggests that phylogenetic effects influence the value of the 
inferred fields h;, and that a correction should be included for 
predictions of energy or fitness. This form of a correction is 
sensible, as phylogenetic effects should make mutations at 
individual residues less frequent, leading to larger inferred fields. 
For closely related strains such as those studied here experimen- 
tally, any systematic inaccuracies in the energy due to phylogeny 
should be similar in magnitude, and thus differences in energy 
should predict relative replicative fitness fairly accurately. This 
would not necessarily be true, however, for sequences separated by 
many mutations. Further theoretical developments may be needed 
to separate out the contributions of phylogeny and intrinsic fitness 
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Figure 3. Classification of HIV-1 NL4-3 Gag mutants as unfit/fit using predicted E values. Graphs show the ability of E classifiers, predicted 
by regularized Ising (panel A) and Potts (panel B) models, to correctly classify HIV-I NL4-3 Gag mutants into unfit (fiC<0.5) and fit (fiOO.5) categories. 
The measured replication capacities of mutants classified as fit or unfit according to their predicted E values were compared with the Mann-Whitney 
test and p values are shown. 
doi:1 0.1 371/journal.pcbl.l 003776.g003 



from the Ising and Potts model landscapes presented here [25] to 
predict the relative fitness of strains that differ by many mutations. 

In addition to phylogeny, other factors such as host-pathogen 
interactions and pure stochastic fluctuations affect the observed 
distribution of sequences, and could complicate fitness predictions. 
In another work [25] we have investigated these issues by carrying 
out stochastic simulations that aim to mimic the way the samples 
were collected and host-pathogen dynamics. In this paper, for the 
pi 7 protein, we found that the fitness and prevalence were not the 
same. However, the rank order of fitness and prevalence were the 
same as long as the strains being compared were not very far apart 
in sequence space. This is largely because of the diversity of 
immune responses due to diverse HLA types in the human 
population. Additionally, the number of virus particles in single 
infected individuals from whom the virus sequences were extracted 
is large, as is the number of patients from whom the virus samples 
were taken. We find that the one and two-point mutational 
probabilities in the sequence databases have converged [11]; i.e., 
these correlations do not change upon removal of some sequences, 
suppressing the effects fluctuations on the inferred model. 

We also note that some caution should be taken in comparing E 
values for sequences belonging to different proteins. The fitness 
predictions of the Ising and Potts models are unchanged by a 
constant shift in energy for all sequences, thus comparisons of 
absolute energy values are not physically meaningful. Differences 
in energy between two sequences in the same protein, however, 
can be unambiguously interpreted as the fitness ratio of those 
sequences. This is the approach we have taken when examining E 
values from sequences with mutations in pi 7 and p24 together: 
rather than comparing the absolute energies, we compare the 
differences in energy between the mutant and the NL4-3 reference 
sequence in each protein, which reflect the fitness of the mutant 
relative to the NL4-3 reference sequence. Finally, translation from 
differences in energy to diflerences in fitness might depend on the 



specific protein that is being considered. While comparisons of 
energy differences and relative fitnesses of pi 7 and p24 mutants 
performed here exhibit no obvious incongruences, further study is 
needed to confirm the generality of fitness predictions across 

proteins. 

In the case of two mutant pairs (186I295E and 186133 IR) that 
were predicted by the models to have very low fitness, partial 
reversions and/ or additional mutations spontaneously arose in 
culture that restored virus viability. However, with the exception 
of one of the spontaneous mutations (2fiOD) observed in 
combination with 186133 IR that modestly decreased the predict- 
ed energy (increased fitness) in the regularized Ising model but not 
the Potts or original Ising models, the models do not predict lower 
energies (increased fitness) for these mutant pairs in combination 
with the additional mutations arising in vitro. Further, three of the 
spontaneous mutations - 186V, 27 IS and 232M - were not 
observed in the MSA and ther(;f()r(- could not be assessed by the 
Potts model. As a possible interpretation of these findings, we 
suggest that it may be the case that these mutation patterns 
observed in vitro are not typically observed in vivo, perhaps since 
these are infrequently explored mutational routes. As a corollary, 
this could indicate an inherent limitation of computational models 
derived from clinical sequence data to identify all possible escape- 
compensatory pathways, and the importance of in vitro and in 
vivo experiments to validate and complement model predictions. 
A mitigating factor, of course, is that mutational pathways 
observed in vitro but not in vivo may be of less direct clinical 
relevance. 

In future \vork, the model predictions will be further validated in 
animal models by testing the viable escape pathways predicted to 
emerge following immunization with immunogens containing 
vulnerable HIV-1 regions only. The validated fitness landscape 
could then be used to design vaccine immunogens containing 
epitopes from the vulnerable regions that could be presented by 
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people with diverse HLAs and that target residues particularly 
harmfial to HIV-l when mutated simultaneously, thereby 
substantially diminishing viral fitness and/or blocking viable 
mutational escape [11,12]. Such immunogens potentially repre- 
sent good therapeutic vaccine candidates to overcome the 
challenge of HIV-1 evasion of CD8+ T cell responses. However, 
further work will also be required to optimize design of such 
immunogens to ensure that epitopes included are processed 
effectively and that they are sufficiently immunogenic, as well as to 
test their immunogenicity, optimal delivery methods and protec- 
tion efficacy in animal models. Furthermore, fitness landscapes of 
HIV-1 proteins may also be more widely applied to identify 
effective antibody targets, and help design potent combinations of 
neutralizing antibodies for passive immunization as well as small 
molecule inhibitors for therapy. 

Supporting Information 

Code SI Source code for the inverse Potts algorithm. A 

compressed file containing source code and instructions for its 
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